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ABSTRACT 

To understand and interpret the observed Spectral Energy Distributions (SEDs) of starbursts, the- 
oretical or semi-empirical SED models are necessary. Yet, while they are well-founded in theory, 
independent verification and calibration of these models, including the exploration of possible degen- 
eracies between their parameters, are rarely made. As a consequence, a robust fitting method that 
leads to unique and reproducible results has been lacking. Here we introduce a novel approach based 
on Bayesian a nalys is to fit the Spitzer-IRS spectra of starbursts using the SED models proposed by 
iGroves et al.l (|2008[ ). We demonstrate its capabilities and verify the agreement between the derived 
best fit parameters and actual physical conditions by modelling the nearby, well-studied, giant Hn 
region 30 Dor in the LMC. The derived physical parameters, such as cluster mass, cluster age, ISM 
pressure and covering fraction of photodissociation regions, are representative of the 30 Dor region. 
The inclusion of the emission lines in the modelling is crucial to break degeneracies. We investigate the 
limitations and uncertainties by modelling sub-regions, which are dominated by single components, 
within 30 Dor. A remarkable result for 30 Doradus in particular is a considerable contribution to its 
mid-infrared spectrum from hot (« 300 K) dust. The demonstrated success of our approach will allow 
us to derive the physical conditions in more distant, spatially unresolved starbursts. 
Subject headings: SED fitting - Hn regions - infrared: ISM - ISM: individual: 30 Dor - stars: 
formation 



1. INTRODUCTION 

In theory, the spectral energy distribution (SED) of 
a galaxy contains a wealth of information about both 
its evolutionary history and current conditions. How- 
ever, extracting this information is difficult and re- 
quires the use of physically based models. Neverthe- 
less, SED fitting is a necessary process as many high 
rcdshift galaxies remain unresolved by our current in- 
struments and any attempts to characterize the condi- 
tions and processes that lead to their starburst activities 
rely almost exclusively on their spatially averaged prop- 
erties. These models of the integrated SEDs of galaxies 
currently cover a wide range of galaxy types, but are 
part icularly dominated by models of Starburst galax- 
ies dGalliano et ail 120031: iSiebenmorgen fc KrugeUl2007t 
i Takagi et alJl2003USilva et al.l 119981; IDopita et aHl2005l 
l2006bild : IGroves et al.ll2008[ ) . The ultraviolet (UV) to far 
infrared (FIR) SED of these Starbursts is dominated by 
the energetic photons emitted by massive stars with typ- 
ical lifetimes of less than 10 Myr. 

In particular, the mid-infrared portion of the SED 
contains several important diagnostics that probe the 
physical conditions of starbursts. Observations of a set of 
marginally resolved starburst galaxies with the Spitzer 
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Space Telescope show a broad range of mid-infrared 
properties, including different strengths of the polycyclic 
aromatic hydrocarbon (PAH) bands, thermal contin- 
uum slopes, depth of the silicate absorption features 
at 10 /im and 18 /zm and intensity of nebular em ission 
lines (jBrandl et al.ll2006t iBernard Saras et al.ll2009l ). All 
these signatures have contributions from different spa- 
tial regions, depending on the geometrical distribution 
of gas an d dust with respect to the ionizing stars. For 
example, Beir ao et al.l (|2009r i reported on the presence 
of compact star forming knots around the nucleus of 
the starburst galaxy Arp 143, and similar star form- 
ing knots have been reported near t he nucleus of NGC 
253 (jFernandez-Ontiveros et al"1 l2009). In other galaxies, 
such as M51, star formation spreads more uniformly over 
the galactic disk. The different distributions of gas, dust, 
and stars in galaxies affect the shape of the spatially inte- 
grated SED. Inversely, a sophisticated and well calibrated 
SED model should be able to recover the information on 
the local starburst conditions from the integrated SED. 

A considerable amount of SED model libraries can be 
found in recent literature (see e.g. iWalcher et al.l [20TTl 
for a comprehensive review on SED fitting). These mod- 
els generally make assumptions on the internal physics 
of galaxies and predict the output SED as a function of 
certain model parameters, such as star formation rates 
(SFRs), metallicity (Z), and the interstellar medium 
(ISM) pressure, density, and temperature, among many 
others. SED fitting refers to the process of choosing from 
a particular library the model solutions that best repro- 
duce the data. While finding the best-fit model via, for 
an example, a \ 2 minimization provides an estimate of 
the parameters, this method alone is insufficient to pro- 
vide absolute parameter uncertainties. In order to obtain 
robust parameter estimates, including uncertainties, it 
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becomes necessary to explore the whole parameter space 
and perform a statistical study of their correlations. We 
highlight four aspects that make this task difficult. First, 
the sensitivity of photometric and spectroscopic studies 
is limited not only by instrumental constraints, but also 
by more fundamental constrains such as shot noise in the 
case of weak sources. Hence, the robustness of SED fit- 
ting depends on the data quality and on sufficient data 
coverage. Second, degeneracies between model param- 
eters are common, especially when limited to a narrow 
spectral window (e.g., the mid- infrared). Third, inde- 
pendent determinations (from observations or theory) of 
the physical parameters against which we can confront 
our model results are rare for most starburst, hence mak- 
ing it difficult to calibrate the models. And last but not 
least, no robust fitting routine that leads to reproducible 
results has been established so far for the specific case of 
starburst spectra. 

In this paper we present a Bayesian fitting routine for 
the mid-infrared (5 — 38 fim) spectra of starbursts that 
can be extended to other wavelengths. We derive prob- 
ability distribution functions (PDFs) for the model pa- 
rameters, and study the implications on the physics of 
starbursts. To calibrate this routine we apply it to the 
mid-infrared spectrum of the 30 Doradus region in the 
Large Magellanic Cloud. The selection of this nearby 
starburst as a calibrator is natural, since its proximity 
(« 53 pc) allows us to differentiate spatially resolved sub- 
regions of the giant Hn region, and study their spectra 
separately. The well studied stellar populations, ionized 
gas, and dust content provide the necessary independent 
measurements to compare with SED fitting results. 

Current spatial resolutions achieved with the mapping 
mode of the Infrared Spectrograph on board the Spitzcr 
Space Telescope are of the order of a few arcseconds at 
5 /xm corresponding to a scale of about one parsec at 
the distance to 30 Doradus. Even the next generation 
spectrometer operating at these wavelengths, the Mid 
Infrared Instrument (MIRI), on board the 6.5 m James 
Webb Space Telescope, will not be able to resolve typical 
giant Hn regions in galaxies located at distances larger 
than about 30 Mpc at a nominal wavelength of 1 5 /xm. 
This highlights the importance of understanding the in- 
tegrated SEDs of these objects. 

This paper is structured as follows. In $2]we describe 
some general aspects of the 30 Doradus region, focus- 
ing on its stellar content and its physical properties, as 
obtained from HST and Spitzer observations, and we dis- 
cuss the Spitzer-IRS spectral data that we model. In $3] 
we give a brief overview of the models we use to generate 
our grid of synthetic SEDs. In S|4] we introduce our fit- 
ting routine and discuss the assumed priors and involved 
uncertainties. $5] presents the results of applying our fit- 
ting routine to 30 Doradus, discuss the implications of 
the model parameters and the physical interpretation of 
the mid- infrared SEDs. Finally, in SjB]we summarize our 
main findings. 

2. THE 30 DORADUS REGION 

Our choice of 30 Doradus as a calibrator relies on three 
powerful reasons: (i) it is the largest giant Hn region in 
the Local Group, (ii) it is well studied across the whole 
electromagnetic spectrum, and (Hi) it is close enough 
to be well resolved into individual components. In this 



section we describe the general properties of 30 Doradus 
and the spectral data that we model. 

2.1. Properties of the 30 Doradus Region 

30 Doradus is the most massive giant Hn region 
in the Local Group. It is located 53 ± 3 kpc away 
(jFeast fc Catchpolelll997| ). in the north-east part of the 
Large Magellanic Cloud (LMC) and includes the stellar 
cluster NGC 2070, the cloud of ionized gas created by 
the ionizing radiation from NGC 2070 and dominated by 
its compact central core R136,and the photon-dissociated 
regions and molecular material associated with the star 
forming region. We show the complexity of the region in 

Fig.[U 

R136 is the most dense concentration of stars in the lo- 
cal group, with an estimated stellar mass of 2 x 1 4 Mq 
contained within the innermost 5pc (jHunter et al.lfl995l) . 
The associated H n region has an Ha luminosity of 
1.5 x 10 40 ergs" 1 (K ennicuttl 119841 ) and a far-infrared 
luminosity of 4 x 10 7 L (|Werner et al.lll978t ). Stellar 
winds, supernovae, and radiation pressure from the cen- 
tral cluster have excavated an expanding ionized bubble 
and created a complex filamentary structure (Fig. [1]). 
This bubble, and other similar cavities in the region 
are filled with X-ray emitting gas at temperatures of 
~ 10 6 K, as revealed by observations with the Chan- 
dra Space Observatory (jTownslev et alJ l2~006). A recent 
study of the optical emission lines shows no evidence 
of ionizatio n by supernova-driven shocks found by a re- 
cent study ([Pellegrini et alJ feOlO). and hence the domi- 
nant excitation mechanism in the 30 Doradus region is 
photoionization by the UV photons produced mainly in 
R136. This was corroborated by a comparison of ob- 
serve d IRS line fluxes with m odels of the mid-infrared 
lines (|Indebetouw et al.ll2009T ) 

Using HST spectroscopy, iWalborn fc Blades! (|1997| ) 
identified several non-coeval stellar populations in the 
30 Doradus region, and classified them as follows: (i) a 
core- ionizing phase (R136), with an age of 2-3 Myr; (ii) 
a peripheral triggered phase, with an age of < 1 Myr 
(this population has also been i dentified using near 
infrar ed excess measurements, e.g. Maerckcr fc Burtonl 
2005); (Hi) a phase of OB supergiants with an age of 
6 Myr; (iv) the Hodge 301 cluster, 3' NW of R136, with 
an age of 10 Myr, and (v) the R143 OB association, 
with ages between 4-7 Myr. 

An interesting aspect of 30 Doradus is its structure 
of bubbles and filaments. Observations of galactic and 
extragalactic H II regions have revealed expanding struc- 
tures of ionized gas driven by stellar winds and supernova 
activity from the OB stellar population. In the particular 
case of 30 Doradus, expanding supershells have been de- 
tected with diameters betwee n 2 and 20pc and expan sion 
velocities of 100-300 km s" 1 (|Chu fc Kennicuttl [l99l . 

The metallicity of 30 Doradus and of the LMC in gen- 
eral is sub-solar (Z = 0.4 Zq) (|Westerhmdlll997D . Due 
to this low metallicity environment, the dust-to- gas ratio 
in the LMC is about 30% l ower than in the Milky Way 
(see review by[Drainc 2003, and references therein), and 
the system allows us to investigate the effect of UV radi- 
ation in lower metallicity environments as compared to 
our own galaxy. 

For simplicity, in this paper we refer to 30 Doradus 
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Fig. 1. — The 30 Doradus region imaged in the 4 Spitzer-IRAC channels. The filamentary structure and bubble-like cavities are evident. 
The ionized gas illuminated by R136 (green) is confined to a thin layer next to the PDR (red), where the PAH emission is found. 

as the region of ~ 100 pc= 4.1 arcmin in diameter in 
projection centered in R136. 

2.2. The Integrated Mid-IR Spectrum of 30 Doradus 

The Spitzer-IRS spectral da ta that we model here ha s 
been extensively discussed in llndebetouw et al.l (2009), 
as part of the Spitzer General Observer Program Stel- 
lar Feedback on Circumcluster Gas and Dust in 30 Do- 
radus, the Nearest Super-Star Cluster, (PID 30653, P. I. 
R. Indebetouw). It consists of four data cubes obtained 
by mapping the 30 Doradus region with the two low- 
resolution slits of the IRS ("short- low" and "long-low") 
in each of their two spectral orders. For reference, the 
first order of the short-low (SL1) map covers an area of 
116 pc x 84 pc, and includes a significant portion of the 
30 Doradus emission nebula. The wavelength coverage 
is between 5-38 /mi with a resolving power R = A/ A A, 
varying from 60 at the short wavelength end to about 
110 at the long wavelength end. Exposure times were of 
the order of 150 s per slit position. 

Spectra of chosen regio ns are extracted us ing the CU- 
BISM software package (|Smith et al.ll2007D . Once the 
sky subtraction has been performed, we extract individ- 



ual spectra using a resolution element of 2 x 2 SL1 pixels 
for all orders. This corresponds to an angular resolu- 
tion of 3.7 arcseconds, and a physical spatial resolution 
of roughly 1 pc at the distance to the LMC. To create the 
spatially integrated spectrum of 30 Doradus, we co-add 
the spectra of all individual resolution elements within 
an area of about 64 pc x 63 pc (the magenta square in 
Fig. [2J . We show the resulting integrated spectrum in 
Fig. El The integrated spectrum is dominated by emis- 
sion from nebular lines and the thermal continuum, while 
the PAH emission is generally weak in the region. 

Here we express all fluxes as vF„ in units of ergs -1 . To 
convert from the MJysr -1 units from the IRS pipeline, 
we multiply the fluxes by the aperture area of 13.7 
arcsec 2 , a nd assume a distance to 30 Doradus (LAIC) 
of 53 kpc (|Feast fc CatchpoldH997l ). 

2.3. Individual Sources 

In Fig. [2] we have indicated four locations defined in 
Table [TJ of which we show their respective spectra in 
Fig. El These locations include sources of different na- 
ture and were chosen to cover a broad range of physical 
conditions and spectral shapes. We model their spectra 
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Fig. 2. — Multi-wavelength view of the 30 Doradus region. Red: 
IRAC 8 fim image showing the PAH emission from the PDR re- 
gion (c.f. [lj. Green: [S IV]10.5 (im emis sion line map, constructed 
from the spectral map described in i|2.2l tracing the distribution of 
highly ionized gas. Blue: Red continuum image showing the stellar 
continuum emission. White circ les mark the positions of the indi- 
vidual spectra discussed in i|2.3l and their sizes correspond to the 
size of one resolution element of the spectral map. The magenta 
square outlines the full IRS spectral map explored in this paper. 
North is up and east is to the left. 
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Fig. 3. — Integrated mid-infrared spectr um of 30 Doradus (red 
line) and spectra of sources described in £|2,3I (black lines) as la- 
belled. All spectra are normalized to the flux at 30 jira and are 
shifted one decade in flux for comparison. The main spectral fea- 
tures are labelled. 



separately to study the validity of the models in environ- 
ments which are dominated by either highly ionized gas 
or by embedded stars. 

Source 1 corresponds to the location of the young OB 
cluster R136. The emission here is dominated by UV and 
optical photons and shows little infrared emission from 
PAHs. 

Source 2 is a YS O candidate selected from IRAC colors 
(iKim et al.l 12007ft , ac cording to the criterion suggested 
by IM en et al about 1 arcminute southwest of 

R136, at the ionized southern edge of the main bubble- 
like structure, in a region with significant [S IV] 10.5 fim 
emission. Its spectrum has a smooth thermal contin- 
uum with no sign of PAH emission, bu with the typi- 
cal nebular lines [Ne H]12.81/zm, [Ne III]15.56^m, and 
[S III] 18.71 /zm. 

Source 3 is a bright infrared source outside of the main 
bubble, to the north-west of the cluster. Its spectrum 
shows prominent PAH emission features and a deep sili- 
cate absorption feature at 10 fim. 

Source 4 is a n infrared source identifie d as a proto- 
stellar object by IWalborn fc Blades! (|1987D . just outside 
the main bubble, north of R136. It coincides with a 
strong peak of [ S IV] e mission and is also an X-ray source. 
lLazendic et al.l (|2003l ) even consider this source to be a 
supernova remnant, but also point to its higher Ha/H/3 
ratio and the possibility of it being an H n region with 
an extinction higher than average. 

In general, we observe that emission from all PAH 
bands is weak towards 30 Doradus as compared to 
other starburst sy stems (see, fo r exam ple the starburst 
SED template in iBrandl et alj (|2006D ). In particular, 
the 17 fim PAH complex generally associated with out- 
of-band bending modes of large neutral PAH grains 
(|Van Kerckhoven et all !2000: Pc eters et al.ll2004f ) is only 
marginally detected in our spectra. A remarkable result 
regarding this point is that the 17 fim complex is weaker 
towards source 3 than expected from the proportional- 
ity relations that have been empirically d erived between 
different PAH bands (| Smith et alJ 120071 ). This propor- 
tionality implies that in starburst galaxies the equivalent 
width of the 11.3 fim feature is about twice the eq uiva- 
lcnt width of the 17 fim feature (|Brandl et al.l 120061 ). If 
this were to hold also for our source 3, we would expect 
a flux density of the 17 fim 20% higher than the thermal 
continuum at this wavelength. However, our data indi- 
cates an upper limit for the 17 fim emission of only 2% 
above the continuum level. 

This suppression of the 17 fim band can have several 
interpretations. A possibility is that the PAH molecules 
are not neutral in this region of 30 Doradus. However, 
source 3 is outside of the main ionized bubble shown in 
Fig. [2j and hence we do not expect a high ionization 
state of the PAHs in this region. Metallicity variations 
could also account f or a change in the relative strength 
of the 17 fim feature (|Smith et al.ll2007l ). but even in very 
low metallicity environments an extremely weak 17 fim 
would also imply a weak 11.3/xm feature, which we do not 
observe. We are left with the explanation of grain size 
effects. As mentioned, emission features between 15 fim 
and 20 fim are associated with larg e PAH grains, typically 
conta ining « 2000 carbon atoms (|Van Kerckhoven et al.l 
2000). Wether the conditions in 30 Doradus are unfavor- 
able for the formation of large PAH grains is the mater 
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TABLE 1 

Localized sub-regions in the 30 Doradus Spectral Map. 



Object RA Dec Remarks 

Source 1 5 h 38 m 42.3 s -69° 06' 03.0" R136 

Source 2 5 h 38 m 49.7 s -69° 06' 42.7" YSO candidate 

Source 3 5 h 38 m 56.5 s -69° 04' 16.9" High extinction (no Mm ss 0.60) 

Source 4 5 h 38 m 48.30 s -69° 04' 41.2" Protostar, [S IV] emission 



of a subsequent paper. 

3. MODELLING THE SEDS OF STARBURSTS 
3.1. Literature on SED modelling 

The simplest Spectral Energy Distribution models con- 
sider a starburst as a single spherical H n region sur- 
rounding a central ionizing cluster, use stellar synthesis 
for the stellar radiation and solve the radiative trans- 
fer for dust and gas in spherical geometry. These semi- 
empirical attempts use observatio ns of specific objects , 
such as star-forming d warf galaxies ([Galliano et al.H2003f) 
or nuclear starbursts (jSiebenmorgen fc Kriigell 120071 ) to 
constrain the model parameters. They are successful in 
reproducing the photometry, and to some extent the IR 
spectra of these objects, but are limited to a narrow range 
of physical conditions (e.g., only two orders of magni- 
tude in dust density). Fully theoret i cal mo dels, such as 
the ones proposed by iTakagi et al.1 (|2003l ) , make simi- 
lar assumptions on geometry, dust properties and stellar 
synthesis, and cover a broader range of physical proper- 
ties to model a larger sample of starburst galaxies, but 
ignore spatial variations of the parameters. 

More sophisticated models consider the starburst as a 
collection of individual Hn regions with different ages 
and environments, whose SEDs add up to produce the 
total galactic SED. In the GRASIL models, for example, 
each of these individual H II regions is assume d to have 
different physical properties ([Silva et al.lll998| ). Unfor- 
tunately, they do not allow for the dynamical evolution 
of the expanding shell-like structures such as the ones 
we have described in fl2.ll In the expanding mass-loss 
bubble scenario, the time-dependent radius and exter- 
nal pressure of the H n region are controlle d by the me- 
chani cal luminosity from the newborn stars ([Castor et al.l 
119751 ). and have a strong influence on the shape of 
the SEP, as they control the gas and dust geometry 
(|Groves et al.ll200l . 

None of the existing starburst models simultaneously 
accounts for both the multiplicity of H II regions in a 
starburst system and their time evolution as individual 
Hll regions evolve as mass-losing bubbles. However, the 
model s described in the ser i es of papers iDopita et al.1 
(12001 .IDopita et all (|2006bl ). IDopita et al.l (|2006d) and 
iGroves et al.l ( 20081) (DfeG models hereafter), represent 
a step forward in our theoretical description of star- 
burst systems, by including these two aspects in a self- 
consistent way. Although these models have been suc- 
cessfully applied to the SEDs of a var iety of objects, 
such as brightest cluster galaxies (BCGs) (jDonahue et al.l 
1201 ID . no systematic study of the model degeneracies 
have been presented. In the remainder of this section 
we briefly describe the underlying physics of the D&G 
models, emphasizing the aspects that are relevant for 



our discussion, and connect this description to the con- 
trolling model parameters. For a detailed description of 
the model, we refer the reader to the Dopita & Groves 
paper series. 

3.2. The physical concept behind the model 

The D&G models compute the SED of a starburst 
galaxy as the sum of the SEDs of individual expanding 
H II regions, averaged over ages younger than 10 Myr. 
By this age, over 95% of the total ionizing photons pro- 
duced during the main seque nce stage of the ma ssive 
stars have been emitted (e.g. IDopita et al.ll2006bT ) and 
the non-ionizing UV flux is rapidly decreasing as the OB 
stars evolve off the main sequence into supernovae. Here 
we will test the applicability of the individual Hll re- 
gions that constitute the building blocks of the models. 
Each individual giant H II region evolves in time as a 
bubble expanding into the surrounding ISM, driven by 
the stellar winds and supernova from the central cluster. 
The dynamical evolution is controll ed by the equations 
of motion of the expanding bubble ([Castor et al.1 119751 ) 
and provides the instantaneous distance of the ionization 
front, dust and molecular gas with respect to the central 
cluster. The time-dependent expansion of this mass-loss 
bubble controls the temperature of the dust and the ion- 
ization state of the gas in the Hll region, altering the 
shape of the SED. 

The stellar synthesis code Star burst99 ([Leitherer et al.l 
Il999t iVazquez fc Leithererll2005l) provides the stellar ra- 
diation field for a population of stars at a given age 
and metallicity. The energy output is normalized to a 
template cluster, whose mass is a free parameter of the 
models and can be scaled to any desired value. The 
stellar mass in the cluster is distributed according to a 
Kroupa IMF with a lower cutoff at 0.1 Mq and an up- 
per cutoff of 120 M ([Kroupal 120021). T he photoioniza- 
tion code MAPPINGSm ([Groves I l200l provides a self 
consistent treatment of both the dust physics (photo- 
electric heating, dust absorption and emissivity proper- 
ties, etc.) and the ISM physics, returning both line and 
continuum emission. The dust surrounding the cluster 
is considered to have contributions from three compo- 
nents: a population of carbonaceous grains with a power 
law size distribution; a population of silicate grains with 
the same size distribution; and a population of polycyclic 
aromatic hydrocarbon (PAH) molecules, whose emission 
is represented by a template based on IRS observations of 
NGC4676 and NGC725 2, both interacting galaxies with 
strong PAH emission ([Groves et alJ [20081 ). Stochastic 
heating is taking into account, and the maximum grain 
size of the distribution is a max = 0.16 /xm. 

The radiative transfer is calculated for two physical 
situations. The first one considers the H II region only 
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and follows the UV photons as they traverse the wind- 
blown bubble, heat the dust and ionize the atomic hy- 
drogen until the boundary of the ionization front. The 
second one assumes a covering photo-dissociation region 
(PDR) around the H II region, with a hydrogen column 
density of logJV(H) = 22.0 (cm" 2 ) (an A v w 1 - 2 at 
solar metallicity). The individual model SEDs are cal- 
culated at a discrete set of ages between and 10 Myrs, 
with a resolution of 0.5 Myrs, and the final integrated 
SED is calculated as the age-averaged energy output of 
the process. 

In S}5] we will use both the integrated models and the 
single H n region models to interpret the observed inte- 
grated spectrum of 30 Doradus. This is equivalent to 
assuming two different approaches for the star formation 
history (SFH) of th region: an instantaneous burst of a 
given age, and a constant SFH over the last 10 Myrs. 30 
Doradus, although dominated by the single star forma- 
tion event that created R136, is neither morphologically 
nor spectroscopically a "single" Hn region, and hence 
both of these simplifying assumptions should be tested 
to encircle the problem. 

3.3. Model Parameters 

The global parameters that represent the general as- 
sumptions of the D&G models and that remain fixed 
by construction are those describing the overall geome- 
try, the stellar IMF, the dust properties, and the PAH 
molecules. In the following we describe the parameters 
that are free to vary in the D&G models. To reduce 
our parameter space and focus our analysis, in our fit- 
ting process we will keep a few of these parameters con- 
stant based on previous knowledge of the region. The 
free parameters are: the starburst metallicity (Z), the 
ISM thermal pressure (P/k), the cluster mass (M c \), the 
compactness (C), the PDR fraction (/pdr), and the mass 
contained in embedded objects (Af em b). 

3.3.1. Metallicity 

We fix the value of this parameter to Z w 0.4 Z@, 
which we consider a good average of several estimates 
using, for example, V LT observations of R R Lyrae star 
and Cepheid variables (jGratton et al.ll2004fl or modelling 
of ch emical abundances in the LMC (jRussell fe Dopital 
1992). Metallicity variations are expected for other ex- 
tragalactic starburst environments, but the well estab- 
lished sub-solar metallicity of the LMC helps reducing 
the parameter space here. 

3.3.2. ISM pressure 

This parameter describes the ambient ISM pressure 
that opposes the expansion of the mass-loss bubble. 
From a comparison between FIR line ratios in a sample of 
star-forming galaxies measured with the Infrared Space 
Obser v atory (ISO) and PDR models by Kaufman et al.l 
(|1999D . lMalhotra et all (|2001[ ) derived thermal pressures 
of the order of 10 5 K cm -3 . The spatial resolution 
achieved by ISO implies that, in most cases, this value 
corresponds to the thermal pressure averaged over the 
entire galaxy. While we acknowledge that Po/k = 
10 5 K cm -3 seems high for the average pressure of the 
LMC, we consider it a reasonable estimate near 30 Do- 
radus, where gas densities have been boosted up by ear- 
lier star formation events. On the high pressure end, 



Po/k is constrained by the pressure of the ionized X- 
ray emitting gas inside the bubble excavated by radia- 
tion pressure near R136, whi ch has been e stimated to be 
of the order of 10 6 K cm" 3 (jWangl I1999Q . We thus fix 
Po/k = 10 5 K cm -3 in our models. 

3.3.3. Cluster mass 

The model SEDs scale in flux according to the total 
stellar mass contained in the star clusters. For an age- 
averaged model, averaged over the last 10 Myr, the scal- 
ing relates to the total mass of stars formed during that 
period of time, and hence the derived mass is interpreted 
as a star formation rate (SFR, in Mq yr _1 ), while for a 
model of a single cluster with a given age (our test case), 
the scaling relates to the cluster mass, M c \. For all cases, 
however, it is assumed that stochastic effects within the 
IMF are limited, and that the stellar population samples 
the full range of stellar masses. 

3.3.4. Compactness parameter 

The D&G models introduce the compactness parame- 
ter, C, resulting from the combination of the ISM ambi- 
ent pressure P and the cluster mass M c \. This dimen- 
sionless parameter characterizes the distribution of the 
ISM with respect to the ionizing stars and is based on 
a constant heating flux input to the stars. Intuitively, it 
describes how close the dust is distributed to the ioniz- 
ing stars as a function of the cluster mass and hence it 
controls the temperature distribution of the dust and the 
far-IR shape of the SED. C is proportional to the time- 
averaged cluster luminosity and inversely proportional to 
the time-averaged square of the swept-up bubble radius. 
As described in D&G, we can define the compactness as: 



!ogC = - log 
5 



M 



log 



P/k 



cm 



K 



(1) 



where M c \ is the cluster mass, P is the ambient ISM 
pressure and k is the Boltzmann constant. The pressure 
parameter P/k relates to the ambient thermal pressure 
(or equivalently, the density) of the surrounding ISM. 

3.3.5. PDR fraction /pdr 

As mentioned above, the D&G models explore two 
cases: a fully exposed H II region (i.e. the ISM ends at the 
ionization front), and an Hn region that is completely 
covered by the PDR in projection, with logiV(H) = 
22 (cm -2 ). In reality, a star forming region will have 
a mix of both PDR emission and direct H II emission, 
which we approximate by the combination of the two 
extreme cases, parametrized by the fraction /pdr: 



F, 



HII+PDR 



fpDRF^ UK + (1 - /pdr)^ 



HII 



(2) 



where jtHii+pdr j g ^ c monochromatic flux arising 
from the star forming region, while F^ UR and F™ 1 
correspond respectively to the fluxes calculated for the 
PDR-fully covered case and the Hn region-only case. 
/pdr = 0.0 implies that there is no PDR material left 
around the ionized region, while /pdr = 1.0 implies a 
fully PDR-covered Hn region. In this fully-covered case, 
the PDR absorbs all of the non-ionizing UV continuum 
and re-radiates it at mid-infrared wavelengths. 
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3.3.6. Contribution from embedded objects 

We expect a considerable contribution from a popu- 
lation of massive protostars to the mid-infrared SED of 
Hll regions and starbursts, due to triggered and ongoing 
star formation. At the early stages of star formation, the 
young objects are in a protostar or Ultra-Compact Hll 
phase, deeply buried in dust envelopes. From an obser- 
vational point of view, and given the age resolution of the 
models, these two types of objects are indistinguishable. 
To account f or them, the model s include a population 
of UCHIIRs (jDopita et al.ll2006afl . In terms of the SED, 
these models add a component of hot dust at around 
25 fim. We parametrize this contribution by scaling it to 
the desired mass, M em b- 

3.4. Attenuation by diffuse dust 

The models include an attenuation factor to account 
for additional absorption of UV light by foreground dif- 
fuse dust. This factor is important in the modelling 
of starburst galaxies, where there is significant diffuse 
material along the line of sight but not associated with 
the star-form ing regions. The a d opted extinction curve 
is derived by iFischera fc Dopital (|2005f) and resembles a 
Calzetti extinction law, which is exponential. The incom- 
ing flux is corrected for extinction as: F = i 7 oe -p * <Jott , 
where p is the column density of dust that gives a certain 
Ay, and cr at t is the dust attenuation cross section Based 
on radio continuum observations, iDickel et al.l (|1994[ ) 
find an extinction of Ay = 1.1 m ag towards the 30 Do - 
radus region. In a recent paper, lHaschke et all ()2011l ) 
find a reddening towards 30 Doradus of E(V-I)= 0.43 
mag, corresponding to a similar extinction. We expect 
individual sources to have higher extinction values within 
30 Doradus, with individual protostars having values of 
Ay up to 4.0 magnitudes. Hence for consistency we use 
here an average value of Ay = 2.0. 

4. FITTING ROUTINE 

We introduce here a Bayesian fitting routine for the 
mid-infrared SED of a starburst, either individual star- 
bursts such as 30 Doradus, or entire starburst galax- 
ies. This routine can be easily extended to include other 
wavelength ranges, and can be used for any observed 
spectrum that is expected to be within the defined pa- 
rameter space. We consider each model parameter as a 
random variable with an associated probability distribu- 
tion function (PDF). Rather than just minimizing the 
X 2 value to find the best fitting model, we solve for the 
probability distribution function of each of the model pa- 
rameters. 

In recent years, Bayesian analysis has been used in 
a number of different fields of astrophysics, where an 
attempt was made to reproduce a limited amount of 
data with multi-parameter models. Some of the ap- 
plications of Bayesian methods in the determinati on of 
best fit parameters include photo metric redshifts (I Wo HI 
2009), observational cosmology (jKilbinger et al.l I2010D 
and dusty tori around Active Galactic Nuclei (AGN) 
(|Asensio Ramos fc Ramos A lmeida 2009). 

4.1. Probability Distribution Functions 

We fit the integrated spectrum of 30 Doradus and the 
individual locations in Table [1] using a grid of the D&G 



models parametrized by the quantities described in JJJ In 
determining the best fit we use ^-minimization, where 
the reduced x 2 is given by 

2 _^ (F l -f(p Q ,\ l )) 2 

Xred 2^ DOF X of ' K ) 

i—i 1 

with the sum performed over all wavelength bins Aj. 
The size of these bins is fixed by the wavelength reso- 
lution of the models. Fi is the measured flux for each 
wavelength bin, / is the model-predicted flux at certain 
wavelength A; for a given set of parameters po, o~i is the 
observational error for Fi, and DOF is the number of de- 
grees of freedom, namely the total number of wavelength 
bins minus the number of free parameters in the model. 
Minimizing \ 2 gives us the best fit values for the param- 
eters, but tells us very little about the uncertainties in 
the model and the parameter degeneracies. However, by 
exploring the \ 2 surface over the range of parameters we 
can explore these degeneracies and the robustness of the 
returned parameters. 

If the errors in the parameters can be described using 
a Gaussian distribution, the Bayes theorem states that 
the probability distribution function (PDF) for a given 
parameter or group of parameters (po) can be recovered 
from the reduced \ 2 distribution: 

P(p )= ^Te- 1 / 2 ^, (4) 

The resulting distribution, called the posterior distri- 
bution for that parameter, is the product of the likelihood 
distribution and a modulating probability distribution 
that includes any available a priori knowledge about the 
parameter, that comes from previous observations, the- 
ory, or the experimental setup. This modulating prob- 
ability distribution is called the prior distribution. We 
refer to the adopted prior distributions as the priors of 
our study. 

4.2. Model Priors 

Initially, we introduce bounded uniform priors for all 
parameters of the D&G model. The bounds introduced 
in our priors are predominantly constrained by theory, 
with some constraints from observations. We use a grid 
of 9 x 10 5 model outputs to cover the broad range of phys- 
ical conditions in starbursts. Table [2] summarizes the re- 
sulting sampling for this study, the parameter ranges and 
their resolutions. 

The range of ages is constrained by the typical main 
sequence lifetime of an early type star. For the compact- 
ness parameter C the limits are related to our knowledge 
of star-forming regions: values below logC = 3.5 would 
imply very diffuse (n « 10 4 /T) ISM or stellar clusters, 
far lower than expected for starburst regions, while val- 
ues exceeding logC = 6.5 would lead to very compact and 
massive clusters. The values of stellar mass in the cluster, 
M c i, and the mass contribution from embedded objects, 
-Memb are selected on a logarithmic scale depending on 
the total mid-infrared flux as measured in the spectra. 
The fraction of PDR material, /pdr 7 ranges from a com- 
pletely PDR-free starburst (/pdr = 0.0) to a situation 
where the Hll region is completely hidden by the PDR 
(/pdr = 1.0). 



TABLE 2 

Values adopted by the model parameters. 



Parameter Range Resolution Remarks 

Age 0-10Myr 0.5 Myr 

logC 3.5-6.5 0.5 

/pdr 0.0 - 1.0 0.05 

Afstars 2 orders of magnitude 0.13 dex Adjusted to total flux density 

M em i, 0.8 orders of magnitude 0.05 dex Adjusted to total flux density 



4.3. Uncertainties and Models Resolution 
4.3.1. Sources of observational error 

There are three types of errors contributing to the total 
uncertainty of the measured flux densities: 

• The absolute flu x calibration . Usin g model stel- 
lar atmospheres, iDecin et al.1 (|2004l ) find that the 
lcr uncertainties on the absolute IRS flux calibra- 
tion are « 20% for the SH and LH modules and 

15% for the SL and LL modules. With regard 
to the modelling this error is similar to an uncer- 
tainty in the distance to the object, and affects 
mainly luminosity-based estimates, such as the de- 
rived SFR or stellar mass. 

• The relative flux calibration. This refers to re- 
sponse variations within the given spectral range, 
often from one resolution element to the next one 
and is the equivalent of a "flat field". From the 
typical differences between spectra of high signal- 
to-noise, taken at two different locations within the 
same slit, we estimated this uncertainty to be about 
5%. With regard to the modelling this error limits 
the weight that can be given to individual spectral 
features, and is thus a fundamental limitation to 
the achievable accuracy. 

• Systematic errors due to the specific observing con- 
ditions. This uncertainty includes observational 
jitter, drifts and source (de-)centering, which may 
lead to a wavelength dependent change in the over- 
all SED slope. It also includes the amount of radi- 
ation that is external to the source of interest but 
picked up by the slit, e.g., from the diffuse inter- 
stellar radiation field or another nearby source. 

Flux calibration of a slit spectrum can assume that 
the source is a point source, and multiply by the 
fraction of the point source outside of the slit, 
which corresponds to a wavelength-dependent "slit 
loss correction factor" . Alternately one can assume 
that the intrinsic distribution of emission is spa- 
tially flat, so the same amount of light is lost from 
the slit as re-enters it from a neighbouring point on 
the sky. CUBISM assumes the latter. Neither ex- 
treme is correct, and it results in a systematic flux 
uncertainty that scales nonlinearly with brightness. 

In addition, many adjacent spectral features, mea- 
sured with the low resolution IRS modules, will 
be blended together. This is most evident for the 
blending of the [Ne II] 12.81 /im line and the 12.7/xm 
PAH feature. Some of these systematic uncertain- 
ties vary in time, location on the slit and wave- 
lengths, and are extremely hard to quantify. Hence, 



we do not attempt to quantify them but we need 
to keep these additional uncertainties in mind. 

Furthermore, many distinct spectral features, such as 
the nebular emission lines, contain information on the 
physical conditions of the ISM which is complementary 
to the information that can be derived from the dust 
continuum. One might thus consider assigning these dis- 
tinct wavelengths a larger weight (i.e., smaller error) in 
the fitting with respect to the more numerous continuum 
bins. However, these line fluxes are difficult to model 
very accurately, and a larger weight combined with some 
mismatch between observed spectrum and model may 
dominate the ^-minimization routine and lead to an in- 
correct local maximum in the PDF. 

One could alos consider a very sophisticated fitting 
procedure where each resolution clement gets its unique 
uncertainty (i.e., weight) assigned, taking all the above 
mentioned error contributions. However, we consider this 
approach practically impossible to reach our main goal, 
namely the provision of a reliable modelling procedure 
that yields reproducible results. From the above argu- 
ments it is evident that any total uncertainty (which will 
be used as weight in the \ 2 fitting) has to be larger than 
the error in the relative flux calibration, but will likely 
be smaller than the absolute flux uncertainty. Our tests 
have shown that a flux uncertainty of 10% for all IRS res- 
olution elements leads to meaningful and robust results. 
Hence, we adopt an error of 10% per IRS resolution ele- 
ment. 

However, the above stated values only hold for the 
Bright Source Limit (BSL) of the IRS, which corresponds 
to a S/N ratio of about 10 {IRS Instrument Handbook). 
For dim sources with S/N < 10 , statistical variations in 
the number of detected photons (i.e., shot noise) domi- 
nate the uncertainties, and our 10% uncertainty estimate 
no longer holds. Other noise sources, such as noise from 
the detector read out and dark current come into play, 
and we need a more conservative error estimate. For all 
spectral resolution elements with S/N < 10 we use the 
RMS variations of the spectrum between adjacent posi- 
tions in the spectral map. For each location we extract 
the spectra of the four nearest resolution elements. From 
these five locations we calculate the average and RMS de- 
viations for each spectral resolution element, and replace 
our standard 10% uncertainty for the BSL by the RMS 
value for that spectral element. 

4.3.2. Data rebinning 

There is a difference between the wavelength bin size of 
the models and the resolution element of the Spitzer-IRS. 
The resolving power of the IRS ranges between 60-110, 
and hence a typical resolution element is of the order of 
0.1 /im. In contrast, for the models we have a wavelength 
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step size that increases logarithmically with wavelength, 
and varies from 0.05 /im at 5/im to l^m at 40 /im, but in- 
cludes local variations to resolve important line features. 
Therefore, we have re-binned the IRS spectral data to 
match the lower spectral resolution of the models, by av- 
eraging the fluxes of the data bins corresponding to the 
same model bin. The uncertainty of the resulting bins, 
on the other hand, are calculated as the square root of 
the quadratic addition of the uncertainties of the origi- 
nal IRS resolution elements. This propagation of error 
with the binning also prevents the uncertainties in the 
long- wavelength bins to dominate the fit. 

5. RESULTS 

In this section we present and discuss the results of 
our fitting routine applied to the integrated spectrum of 
30 Doradus, and subsequently also to selected subregions 
within 30 Dor. The latter have been added (see also sec- 
tions 2.3 and 5.3) to probe the validity and limitations of 
our starburst models on regions whose spectra are dom- 
inated by one type of source (e.g., an OB cluster or a 
protostar). 

In the following subsection we compare the results on 
the integrated starburst spectrum of 30 Dor, using three 
different approaches for comparison: (i) fitting all res- 
olution elements of the entire 5 — 38/im spectrum, (ii) 
fitting the continuum bins only, i.e., excluding the fine- 
structure emission lines, and (Hi) deriving the stellar ages 
from the IRS high resolution lines only. We will show 
how sensitive the results depend on the spectral infor- 
mation provided. We then compare the best fit results 
on the integrated 30 Dor spectrum with the results on 
the subregions. Finally, we discuss how the results would 
change if we would not assume a single age burst but an 
age averaged model. 

5.1. Nebular lines ratios as age diagnostics. 

The collection of emission lines in the mid-infrared 
wavelength range of the IRS are sensitive diagnostics 
of ages of massive OB stars and to the hardness of 
the radiation field. They have the advantage of suf- 
fering little from obscuration and hence allow us to 
probe the conditions of deeply buried regions. In 
particular, the [Ne III] 15. 5/im/ [Ne II]12.8/im and the 
[S IV]10.5/xm/[S III]18.7/im line ratios are good diagnos- 
tics of the ionization state of the gas, and hence provide a 
good constraint on the age of the ionizing cluster through 
measuring the har dness of the radiation field (see e.g. 
IGroves et al1l2008h . 

By using the ratios of the nebular emission lines, we 
can provide useful constraints on the ionization state and 
age of the central cluster, and thus break the present 
degeneracies. In the particular case of 30 Doradus, we 
have measurements of the line fluxes with both the low 
resolution orders {lores hereafter) and the high resolution 
orders {hires hereafter) of the IRS. The lores lines have 
a larger flux uncertainty, and hence a good consistency 
check is to compare the results we get from the routine 
with results derived from the hires lines. 

We use the hi r es ne bular line fluxes presented in 
iLebouteiller et all (|2008l) as a reference for the estima- 
tion of cluster age and ionization parameter and compare 
the results with what we obtain from the SED fitting for 
the four individual sources of Table [TJ These line ratios 



ar e plotted in Fig. E] supe rimposed on to a grid of models 
bv lLevesaue et al.l (12010T) (cre ated using the ITERA pro- 
gram of IGroves fe A llen 2010). These are essentially the 
same as the D&G models, and use both the Starburst99 
and MAPPINGSiii co des with similar a ssum ptions about 
the gas. However, the iLevesque et all (|2010l ) models use 
a much simpler geometry (namely plane-parallel instead 
of spherical), demonstrating much more clearly how the 
degeneracy between the hardness of the radiation field 
(i.e. stellar cluster age) and the ionization parameter, Q 
(the ratio of the ionizing photon density to gas density), 
is broken using four strong mid-infrared emission lines. 
The dependence of the line ratios on these two p aram- 
eters has also been note by Morisset et al. (2004). Our 
comparison between the measured line ratios and the pre- 
dictions from the Levesque models indicate ages between 

2.0-2.5 Myr for all four pos itions. 

The sample of sources in ILebouteiller et all (|2008[ ) in- 
cludes five more locations in the 30 Doradus region, apart 
from our four selected sources. Assuming that the mea- 
sured line fluxes in these sources are representative of 
the overall conditions in the cluster, we estimate the line 
ratios for the whole region from the luminosity-weighted 
average line fluxes: log 10 [Nc III]/[Ne H] 3 oDor = 0.75 and 
log 10 [S IV] /[S in] 30 Dor = 0.005. We also plot this average 
value in Fig. @] The result indicates an age of 2.5 Myr. 

I I I 1 2.50 
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Fig. 4. — Measured mid-infrared nebular line ratios ov erplottcd 
on a grid of starburst models from Levesque et al. (2010) that in- 
clude a single H II region only. The parameters of the model are 
the age of the stellar cluster, and the ionization parameter, Q (the 
ratio of ionizing photon density to gas density), as labelled on the 
color-bars to the right. 

In Table|3]we list the ages derived from the Lebouteiller 
et al. hires line ratios as compared to the ages derived 
from the SED fitting in two cases: (a) fitting the emission 
lines and (b) excluding the emission lines. The reason 
to perform the fit using the continuum only is two-fold. 
First, we want to check the consistency of the results 
for our individual sources, which can not be treated as 
isolated H II regions, since their ionization states are af- 
fected by other external sources. Second, as we have 
already stated, while the line ratios can be reasonably 
estimated by the models, the lores equivalent widths of 
the lines are predicted with a lower degree of accuracy 
by the D&G models. 
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TABLE 3 
Ages of individual sources 



Source SED fit, with lines (Myrs) SED continuum (Myrs) Line ratios (Myrs) 



1 7.5 2.0 2.0-2.5 

2 10.0 2.5 1.5-2.0 

3 8.0 2.5 2.0-2.5 

4 0.5 2.0 2.0-2.5 
30 Dor 3.0 5.0 2.5 



The results show that for the individual sources, the 
ages derived from the continuum-only fit are consistent 
with the high resolution measurements from Fig. 0] and 
with the independent measurements of the overall age of 
the region. On the other hand, including the unresolved 
lines in the fit for these individual sources leads to age 
estimates which are in disagreement with all the other 
methods, with a tendency to overestimate the ages. 

The lowres line ratios imply age estimates that are not 
significantly different from the hires results and hence the 
mentioned disagreement should be interpreted in terms 
of the limitations of the integrated H n region D&G mod- 
els to reproduce the continuum and the line emission for 
individual sources. Nonetheless, for the integrated spec- 
trum, the derived age from the emission line fit is con- 
sistent with the high resolution measurements and with 
the literature, as expected for a self-contained region for 
objects of which class the models were intended. 

Based on this results, for the integrated spectrum of 
30 Dor we run the fitting routine including the emission 
lines. For the individual sources, however, we do not at- 
tempt to fit the low resolution lines and fit only the con- 
tinuum. Also for the individual sources, to include the 
information contained in the high resolution line mea- 
surements, we modify the prior probability distribution 
for the ages from those listed in Table [2] to use a Gaus- 
sian distribution centered a 2.5 Myrs with a dispersion 
of 1.5 Myrs. This suppresses weights solution with older 
ages down, further constraining the parameters. 

5.2. Integrated Spectrum 
5.2.1. Best fit 

We show the resulting best fit from our code to the 
integrated spectrum of 30 Doradus in Fig. [5j with the 
residuals of the model fit to the observations shown in the 
lower panel. The quality of the fit is remarkable, with 
most of the spectral features in the mid-infrared spec- 
tral range successfully reproduced. This is a significant 
improvement from broad band photometry SED fitting, 
where only a few data points were fitted to constrain an 
equal number of parameters. 

The individual contributions from the unobscured H II 
region, the PDR, and embedded populations are explic- 
itly plotted in Fig. [5j The residuals in the bottom panel 
indicate that the model fits the observations within the 
uncertainties for most of the IRS wavelength range, but 
underestimate the fluxes near 15/^m. This feature is most 
likely due to an overabundance of small silicate grains 
within the assumed dust model, and is dominated by the 
embedded star model. 

The sulphur lines at 10.5 fim and 18.3 /Ltm appear as 
underestimated by our best fit model. This could be pos- 
sibly due to abundance and/or pressure variations (see 



iDopita et al.ll2006cl ). and we can only argue here that 
for the assumed abundances and ISM pressure, the fit in 
Fig. [5] represents the best case in the prediction of line 
ratios. The mid-infrared continuum is dominated by the 
embedded population, especially for A > 10 /im with the 
PDR contributing mainly to the PAH fluxes and the con- 
tinuum slope at the long wavelength end of the spectral 
range. The Hn region and PDR are responsible for most 
of the emission lines. 
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Fig. 5. — Best fit to the integrated IRS spectrum of 30 Doradus. 
Red: Observed spectrum rebinned to the model resolution with 
error bars for each bin. Black: best fit SED. Dashed blue: "Naked" 
Hn region contribution. Solid blue: PDR and obscured Hn region 
contribution. Dotted blue: Embedded object contribution. The 
best fit values and reduced x 2 are indicated. Residuals are shown 
in the lower panel, in the same logarithmic units. 



5.2.2. Interpretation of the results for the integrated 
spectrum 

The normalized PDFs for the model parameters are 
shown in Fig. |6l for the priors as listed in Table [2] (dot- 
ted lines) , and for the modified probability distribution of 
ages described earlier (dashed lines). The best fit values 
marked by vertical lines and the la uncertainties indi- 
cated by the horizontal line pattern. We list the best-fit 
values, with the uncertainties corresponding to each case, 
in Table H 

From the dark-shaded and dotted PDFs in Fig. [5] it is 
evident that several of the parameters appear to be very 
broad or even unconstrained. The reason for this can be 
clearly seen when we plot 2D PDFs for selected pairs of 
parameters in Fig. [7] (i.e. collapsing the \ 2 space down to 
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two parameters). These show degeneracies between cer- 
tain model parameters, indicating that, at least in the 
IRS wavelength range, these parameters affect the SED 
shape in a similar way. If one or both of these param- 
eters can be constrained using other information, such 
as from other wavelengths, the ID PDFs should become 
narrower, and the parameters better constrained. 

In order to understand these degeneracies we need to 
look carefully at the 2D probability maps and link the 
resulting distributions to the effect that each parame- 
ter has on the spectrum. There is an age-compactness- 
cluster mass degeneracy revealed by two different set of 
parameters that provide a good fit to the observed SED. 
The first panel of Fig. [7] clearly shows the resulting two- 
peak distribution on the probability distribution for the 
total cluster mass-age subspace. These two parameters, 
as well as compactness, have a similar effect on the mid- 
infrared continuum as they vary across the grid: they 
scale the continuum flux by certain multiplicative factor. 

The age-mass component of this degeneracy is not sur- 
prising: as compared to a cluster of certain mass and age, 
a less massive cluster is dimmer at Spitzer wavelengths, 
but the same holds true for an older cluster. IRS contin- 
uum fitting alone is incapable of distinguishing between 
these two parameters, as can be seen from the two-peak 
distribution. However, we have more information con- 
tained in the nebular lines. In particular, older clus- 
ters show less nebular emission, as the ionizing radiation 
strongly decreases with age. Including the lines in the fit 
breaks the age-mass degeneracy and enables us to select, 
between the two possible solutions, the one that best re- 
produces the measured line ratios. The best fit in Fig. [5] 
corresponds to this best solution. 

The two peaks of the compactness-age degeneracy are 
clearly seen in the PDF maps, as shown in the second 
panel of Fig. [7] This degeneracy implies that both a 
young cluster with small compactness or an old cluster 
with high compactness lead to similar fits of the observed 
spectra, provided that the cluster mass also adjusts. This 
is shown in the central panel of Fig. [7] In this case, the 
emission lines are less sensitive to the variations of the 
two parameters together, since both the age of the clus- 
ter and the compactness affect the line ratios. However, 
once the age has been determined from the line ratios, 
the compactness probability distribution also shrinks and 
selects only one of the possible solutions for compactness. 

Panel (c) of Fig. [7] shows a degeneracy between cluster 
mass and PDR fraction. The two strips correspond to 
the two peaks of the age-mass-compactness degeneracy, 
while the smooth diagonal variation corresponds to the 
PDR fraction-cluster mass degeneracy. This degeneracy 
arises from the fact that the PDR region that covers the 
Hn region contributes mostly PAH emission, but also 
adds thermal dust continuum that in the models scales 
up with the PDR fraction. The emission lines are not of 
great help in breaking this degeneracy, since their relative 
fluxes are almost insensitive to variations in PDR content 
and total mass. In this particular point, thus, we can 
only do better if we include data from other wavelengths. 

To assess the importance of the line ratios in the con- 
straining of the parameters and the break of the age- 
compactness-mass degeneracy, in Fig.[5]we plot the PDFs 
for the same parameter pairs, but this time after only the 
continuum has been fitted. A quick comparison between 



the two cases reveals that the inclusion of the lines not 
only selects one of the two degenerate peaks, but also 
helps the best fit values to converge towards the abso- 
lute maximum of the PDF. This is particularly evident 
for the mass-age degeneracy, where the continuum-only 
fit favors an older age solution, while the inclusion of the 
lines shifts the probability maximum to an age that is in 
agreement with independently measured values. 

5.3. Individual sources 

As noted before, any complex starbursting system is 
likely to cover a wide range of object types and physi- 
cal conditions, from individual protostars and luminous 
UCHIIRs to OB clusters, loose stellar associations, PDRs 
and the diffuse ISM. To investigate the range of condi- 
tions for which our starburst models still yield accurate 
results, we have chosen four subregions that probe these 
"extreme" cases where one of these components is ex- 
pected to dominate the mid-infrared spectrum, based on 
detection of, for example, infrared excess or X-ray emis- 
sion. These sources are the OB cluster R136 with little 
dust obscuration, an H II region which shows high ex- 
tinction along the line of sight (source 3), and two com- 
pact objects, which are luminous protostellar candidates 
(sources 2 and 4). 

We fit the continuum spectra of those subregions and 
include the emission line ratios as a modified probability 
of the ages, as described in 35.11 While we do not expect 
to get a very good fit on these types of sources with the 
general starburst models, we want to verify that the cru- 
cial parameters are qualitatively still constrained within 
reasonable limits, according to the respective physical 
condition probed by the individual sources. 

5.3.1. Best fits to the spectra of individual sources 

We have seen that in the case of the integrated spec- 
trum, fitting the emission lines along with the continuum 
greatly helps in breaking the model degeneracies. How- 
ever, as discussed in §5. H and summarized in TableEl this 
is not the case for the individual sources, where the fit- 
ting of unresolved lines leads to age estimates which are 
in disagreement with the continuum-only fit and with the 
line ratio analysis, for reasons that are described in 35.11 
Hence, we do not include the emission lines in the SED 
fitting of the individual sources. Furthermore, we mod- 
ify the age prior to include only ages that are consistent 
with the high-resolution line ratios, within the uncer- 
tainty limits set by the comparison of the line ratios and 
the Levesque models. In each case, instead of the uni- 
form prior distribution of probability for the age, we use 
a Gaussian PDF centered at 2.0 Myr with an uncertainty 
of 1.5 Myr. 

The resulting best-fit SEDs are shown in Fig. HI The 
resulting best fit parameters and \-a ranges derived from 
the PDFs are shown in Table [5] 

5.3.2. Interpretation of the results for individual sources 

All of our spectra show significant flux densities in the 
10//m range. This continuum emission is indicative of hot 
dust at T » 300 K, which is typically associated with pro- 
tostars, but not exclusively. It may also include emission 
from dust close to slightly more evolved stars, as well as 
hot dust in between stars of a cluster, and dense clumps 
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Fig. 6. — Probability distribution functions of the model parameters when fitted to the integrated spectrum of 30 Doradus, for two cases. 
Dotted line: uniform priors described in Table [2] Dashed line: Modified prior distribution for age, as defined in £|5.1I The shaded areas 
correspond to the 1-<t integrated probabilities, while the vertical lines indicate the best-fit values. 



TABLE 4 

Best fit to the integrated spectrum of 30 Doradus 



Parameter Age-unconstrained Age-constrained Literature 
Age (Myr) S.O+^g ~ 3 (Hunter et al. 1995) 

logC b.OtH 5.0l°J 



/pdr 0.6t° : J 0.61°;* 

8 +0.7 , s+0.2 



logMd (Mq) 4 - 8 -o'.i 4 - 8 -a2 4.7 for R136 only (Hunter et al. 1995) 

0.08 
0.02 



l°gM emb (Mq) 4.471°;°® 4.47+ ' 08 



TABLE 5 

Best fit parameters for individual positions in 30 Dor 



Source 1 Source 2 Source 3 Source 4 



t (Myr) 


2 0+ 20 


9 c+2.0 


9 c+2.5 


2.0 


logC 


q c+1.0 


q cr + 0.5 


q tr+0.5 


4.5 


/pdr 


o oo+°- 40 


10+ 050 


25+ - 70 

u - zo -o.oo 


0.80 


logAf cl (M ) 


9 i +0.0 

■^-l.l 


i k+0.0 
1 -°-0.6 


9 9+0.2 
z - z -0.5 


2.1 


logM cmb (Mg) 


73+ 12 


i ,,+0.07 
1 - lo -0.03 


61+ 014 


1.43 



in the H II region that cannot be modelled by the simple 
uniform H II region model of D&G. We account for all 
these contributions by what we have called the embedded 
component. Fig. [5]shows that this component dominates 
the emission of 30 Doradus at mid-infrared wavelengths. 

In fact, our attempts to fit the integrated spectrum 
of 30 Doradus without including this embedded compo- 
nent have proven unsuccessful, and hence, it is one of our 
main results that this component of "embedded objects" 



is necessary to fit the observed spectra for A > 10 ^m. 
Qualitatively, this interpretation looks quite plausible: 
the two positions that coincide with the location of YSO 
candidates (sources 2 and 4) have the higher relative 
mass contribution from the embedded component, with 
30% and 20% of the total mass contained in embedded 
objects, respectively. The corresponding contributions 
from embedded mass in R136 and the highly extincted 
source are 5% and 2%, respectively (Table [5]). However, 
we need to keep two issues in mind: 

First, while we associate this component with a re- 
cently formed star, strictly speaking it is not entire due 
to protostars and UCHIRs, for the reasons given above. 
Hot dust in starbursts may also be found in other envi- 
ronments and, hence, the amount of "embedded objects" 
derived from our fits can only be considered as an upper 
limit on the amount of protostars and UCHIRs. 

Second, the contribution of this embedded component 
to the total emission at A > 10 /xm may be surprisingly 
high but is not unreasonable. This is illustrated in Figs. 
[S]Ja) and[5Jb), which show the fits to source 1 (R136, the 
main cluster) and source 2 (a protostar, or group of pro- 
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Fig. 7. — Two-dimensional PDFs for selected pairs of parameters showing the model degeneracies. The color code indicates normalized 
probability. The cross symbols mark the best-fit values while the color contours indicate the 1-cr (blue) and 90% (green) confidence levels. 
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Fig. 8. — Two-dimensional PDFs for selected pairs of parameters, when only the continuum has been fitted. The color code indicates 
normalized probability. The cross symbols mark the best-fit values while the color contours indicate the l-cr (blue) and 90% (green) 
confidence levels. 



tostars) . While M c \ (stars + H n region + PDR ensem- 
ble) is at least one order of magnitude larger in the case 
of R136 than in the embedded region, the contribution 
of the embedded component is much higher for the pro- 
tostar in comparison, and so are the observed flux den- 
sities over most of the IRS spectral range. While R136 
contributes most of the stellar mass, much less massive 
components described as "embedded objects" contribute 
the majority of the mid-IR flux. 

In other words, the fact that the integrated flux in 30 
Doradus is dominated by this embedded component does 
not imply that there is a similar mass contribution from 
this embedded component. In fact, our derived mass 
contribution of embedded objects to the stellar mass of 
30 Doradus is about 35%. Furthermore, despite the un- 
certainty in the nature of the embedded component, we 
will show in section 5.4 that the star formation rate in 30 
Dor derived from our modelling approach does not over- 
estimate the "true" star formation rate as derived from 
integrated panchromatic SEDs that include the far-IR. 
Additionally, the application of our routine to starburst 
galaxies shows that this component of embedded objects 
does not dominate the mid-infrared emission for these 
galaxies, as it does for 30 Doradus. 

Our results in Table [5] indicate that source 4 has a 
higher compactness as compared to the other individual 
sources. High values of logC are expected in compact 
starbursts with high surface brightness, where dust is 
in close proximity to intense UV fields. Source 4 is a 
very bright and compact source of [S IV] 10.5 /iin, indi- 
cating the presence of highly ionized gas probably near 
a hard UV source, and, as pointed out in ^2.31 is also 



a bright X-ray source that has even been considered as 
a supernova remnant candidate. No other location in 
our spectral map shares these characteristics as an indi- 
vidual source. On the other hand, low compactness is 
derived for sources where the simultaneous presence of 
bright stars and dust can be inferred, as it is the case of 
sources 1, 2 and 3. 

We conclude that, even though our models are not 
intended to model these individual sources, we can 
nonetheless learn from the compactness parameter, as 
defined in §3.21 by comparing the results of the routine 
applied to them. The routine is capable of constrain- 
ing the proximity of luminous sources and hot dust. We 
have derived a relatively high compactness for 30 Do- 
radus itself with logC = 5.0, which is consistent with 
the luminous cluster in its center surrounded by nearby 
ridges of dust. 

A comparison of the values for /pdr in Table [5] with 
the spectra in Fig. [3] indicates that, as expected from 
the construction of the models, a high covering fraction 
is generally associated with strong PAH features. We 
show here that source 2, where a YSO has been identi- 
fied via infrared excess, has little associated PAH emis- 
sion, whereas source 4, which also shows infrared excess 
and has been associated with a YSO, shows significant 
PAH emission. The presence of PAH emission in the 
line of sight towards embedded objects might depend on 
the evolutionary stage of the YSO, i.e., on the optical 
thickness of the envelope. If UV photons from a young, 
massive protostar manage to escape the embedded region 
and create a PDR around the YSO, then we expect to 
detect PAH emission. In some cases, like source 3, where 



14 





Fig. 9. — Best fit models for the individual positions, excludingthe emission lines from the fit and with the age distribution constrained 
by the line ratio analysis. The color code is the same as for Fig. [5] 



the deep silicate feature typical from highly embedded 
objects is accompanied by high line ratios indicative of 
an ionizing source, the PAH emission can even domi- 
nate over the embedded component. No significant PDR 
emission is inferred from the fit to R136. This is con- 
sistent with the expected absence or low abundance of 
PDR material very close to the ionizing cluster. 

There are two important caveats that we must consider 
in the interpretation of the PDR covering fraction. First, 
we have not included here emission from diffuse dust, not 
associated with the starburst (i.e. heated by the inter- 
stellar radiation field); this adds both cold dust and PAH 
emission. Second, there exists a degeneracy between the 
covering fraction of PDR material and the cluster mass 
in the mid-infrared (right panel of Fig. [7]) , due to the fact 
that PDR regions not only add PAH emission, but also 
continuum emission at longer wavelengths. This trans- 
lates in the apparent mismatch between the PDFs and 
the best fit values for these two parameters in Table [5l 

5.4. Age averaged case 

The degeneracy found for the integrated spectrum of 
30 Doradus that leads to the two possible solutions of 



an old, massive cluster or a young cluster with a smaller 
mass is indicative that not even 30 Doradus can be con- 
sidered as a single coeval stellar population. As described 
in the introduction, a number of spectroscopically identi- 
fied populations have been identified in 30 Doradus, and 
the fact that the continuum is compatible with two dif- 
ferent sets of parameters leads us to the conclusion that 
a simple, single-age approach might be inaccurate even 
for our benchmark Hn region. Thus we also carry out 
the SED fitting using an age-averaged model. 

We have fitted the integrated spectrum of 30 Dor us- 
ing the age-average model described in M3.21 where the 
SED is integrated over the 0-10 Myr lifetime of the ion- 
izing stars. In this case, the cluster age is no longer a 
free model parameter, and we can interpret the absolute 
flux scaling as the SFR instead of a single cluster mass. 
We keep the ISM pressure (P/k) and the metallicity (Z) 
fixed to the same values as for the single age case. The 
best fit parameters we calculate in this way must be in- 
terpreted as average values over the time span covered 
by the models. Fig. [10] shows the resulting best fit to 
the data with the age-average model, and Table |6] lists 
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TABLE 6 

Best fit to the integrated 
spectrum of 30 doradus with 
age-average model 



Parameter 



Best fit 



logC 
/PDR 

logSFR cff (M yr- 1 ) 

/cmb 



r o+O.O 

°- u -i.o 
i oo+ - 00 

z - uo -0.05 

o on - 64 



the best-fit parameters with the associated uncertainties 
derived from the PDFs. 




Metallicity = 0.4 Z 
Log C = 5.0 
W = 1-00 
SFR elf = 0.009 M yr" 

femb = 3.16 

P/k = 5.00 



0.4 1- I 



-+- 



10 15 20 25 30 
Wavelength (//m) 



35 



Fig. 10. — Best fit to the integrated IRS spectrum of 30 Doradus 
using the age-average model. The color code is the same as in 
Fig.fi] 



Here we use the parameter / cm b instead of M em b- It 
refers to the ratio of mass contained in embedded objects 
to the mass of main sequence stars for objects younger 
than 10 Myr, and hence it is related to the amount 
of currently ongoing star formation. In other words, 
/cmb gives the fraction of embedded/UCHII luminosity- 
weighted contribution that we have to add to the SED 
to fit the observed spectrum. If / cm b = 0.0, there is no 
current star formation happening, whereas if / cm b = 1.0, 
half of the massive stars formed over the last million year 
are still in a embedded state. Since this contribution is 
integrated over a period of 1 Myr only, adding embed- 
ded objects also implies that the average SFR has to be 
modified according to: 



SFR eff = SFR + ^ SFR 

IAt 



(5) 



where SFR c ff is the effective SFR that accounts for 
the additional population of embedded objects, and /a* 
is the ratio of the total time over which the starburst 
has been modelled (10 Myr) and the estimated duration 
of the embedded phase (1 Myr). For our best fit case, 
we get SFR = 0.007 M yr" 1 and / Gmb = 3.2, which 
implies that among the stars younger than 1 Myr, there 



are about three times more embedded objects than main 
sequence objects. Hence, the effective SFR is SFR c g = 
SFR + 0.32 x SFR = 0.009 M yr" 1 . 

Th e LMC has a SFR of 0.1 M yr" 1 (|Whitnev et all 
2008). Keeping in mind that this is only a lower limit es- 
timate, given the incompleteness of any YSO catalogue, 
our result implies that between 3% and 10% of the star 
forming activity of the LMC takes place in the 30 Do- 
radus region. To see how this compares to estimates of 
the SFR in 30 Doradus from single photometric measure- 
ments, we compare the IRAS flux at 25 y«m for the entire 
LMC to the 64 pcx63 pc area from which we have ex- 
tracted the spectrum of 30 Dor. The total flux density 
from the LMC at 25 /mi is 7520 ± 1100 Jy ([Israel et all 
2010). From our integrated spectrum (Fig.[3]), we derive 
a monochromatic flux density of 1739 ± 174 Jy for the 
same wavelength, which corresponds to 24% of the to- 
tal LMC flux density. Assuming this wavelength directly 
traces star formation, it suggests that our value is close, 
but may be underestimated by a factor of two. 

The compactness and fraction of PDR results are con- 
sistent in both the single age case and the age average 



case. 



6. SUMMARY AND CONCLUSIONS 



Significant progress in our understanding of starburst 
systems has been made over the past decades, both on 
the observational and theoretical side. A huge amount 
of spectral data on star forming regions and starburst 
galaxies has been collected with ISO, Spitzer and the 
Herschel Space Telescopes, complemented by a consider- 
able library of SED models that predict the energy out- 
put of starbursts as a function of wavelength. 

However, SED fitting of starburst has mainly focused 
on maximum likelihood methods, which generally over- 
look degeneracies between physical parameters and lead 
to results that are not unique. Furthermore, these ad 
hoc approaches often depend on some hidden assump- 
tions that make the results reproducible. In this pa- 
per, we presented a routine to fit the SEDs of starbursts 
based on the mod el s proposed in the se r ies of papers 
iDopita et all (120051) . D opita et all (|2006b| ). lDc~pita et al.l 
(|2006cD and iGroves et all (|2008lL We verified the ac- 
curacy and limitations of our approach by comparison 
between the model fit results and the known properties 
of the well-studied, prototypical giant H II region 30 Do- 
radus. Our main findings are: 

• Our modelling procedure is able to fit a broad range 
of continuum slopes, PAH intensities, and emission 
lines. Although we have only used the mid-infrared 
spectra for the calibration, the method can be eas- 
ily expanded to other wavelength ranges. 

• We have verified the validity of our approach by 
comparison with the well studied 30 Doradus re- 
gion. The derived physical parameters, such as 
cluster mass, cluster age, ISM pressure and PDR 
content, are in good agreement with the known 
properties of this nearby starburst. 

• We have provided a detailed study of the model de- 
generacies in the mid-infrared window of the spec- 
trum, and have shown that the best fit values to the 
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continuum shape are driven by a triple luminosity- 
age-compactness degeneracy that, in general, leads 
to multiple "best fits" . 

• The inclusion of emission lines in the analysis 
breaks this degeneracy. It is expected that the 
addition of other wavelength ranges would further 
constrain the model parameters. In particular, the 
precise location of the peak of the dust emission 
in the far-infrared is crucial to constrain the com- 
pactness parameter. Herschel spectroscopy, as well 
as MIPS and PACS photometry, play an important 
role here. 



■ We provided meaningful results to the model- 
defined compactnes s parameter C, introduced by 
IGroves et al.1 (|2008f ). and linked them to the prox- 
imity of ionizing sources and hot dust. 



• We have shown that modelling the SED of a typical 
starburst region requires a component of heavily 
embedded objects (massive YSOs and UCHIIRs) 
which dominate the mid-infrared continuum slope. 
The derived mass fraction of this embedded com- 
ponent can be interpreted as an upper limit to the 
amount of current star formation, since there are 
other dust heating mechanisms not included in the 
models. 



• We found a degeneracy between the total stellar 
mass and the relative amount of PDR material, 
/pdr, as both will contribute to the dust contin- 
uum. This degeneracy may lead to uncertain mass 
estimates and can only be resolved with additional 
data at longer wavelengths, e.g. from Herschel. 

• Generally, two critical assumptions in all starburst 
models are the age and duration of the burst. 
Our "local" template 30 Doradus nicely illustrates 
the typical complexity of a starburst with both, 
the presence of a luminous, coeval cluster (R136), 
and strong evidence for continuous star formation 
across the region. Hence, we have also used an 
age-average model of continuous star formation for 
comparison. This model delivers values for com- 
pactness and PDR contribution that are consistent 
with those derived from the single age models. For 
30 Doradus we derive a contribution of approxi- 
mately 10% to the total SFR of the LMC. 

Now with a robust and well-tested modelling and fit- 
ting routine in hand, we plan to apply this approach to 
more distant giant Hn regions and starburst galaxies. 
The lack of spatially resolved data on e.g., more distant 
ULIRGs and sub- millimeter galaxies requires reliable and 
well calibrated models to derive the physical conditions in 
these starbursts. The novel fitting procedure presented 
in this paper constitute the next step in starburst mod- 
elling and puts such studies on solid grounds. 
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